The following three m-files are a little more elegant than the ones in the "My Big Saab Jump" article.
As provided, they will run on Octave version 4.2.1 and FreeMat version 4.0. They can be easily modified
for other vehicles or other projectiles. 

To use these files, copy all of the text between dashed lines into your favorite text editor and save each
to a directory of your choice using the name provided in the first dashed line. The first file, the main script,
can be saved under any name of your choosing as long as it retains the m extension. However, if you change the
name of either of the other two files, the first script must be changed to reference the new name.

File #1, main script
----------------------------------------------- MySaab96Jump.m -----------------------------------------------
%% MySaab96Jump.m, Ver. 5, 3/31/2018
%%
%% Copyright 2018, George (Ted) Yurkon
%% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. To
%% view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative Commons,
%% PO Box 1866, Mountain View, CA 94042, USA.
%%
%% In short, this work may be shared provided terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%% International License are adhered to.
%%
% Based on projectile theory by Colin Ratcliffe at...
% http://ratcliffe.site/USNA/EM375/labs/07Project/ProjectileTheory.pdf
% Assumed units: meters, seconds, Newtons, kg, radians
% This program runs successfully on FreeMat v4.0 and Octave 4.2.1
% It may run on Matlab as well, but will need modifications
% for Scilab. Scilab provides a Matlab to Scilab converter which may help.
%%
% Some enhancements to the plot figure are available but Octave and 
% FreeMat are not compatibe. Seach for "---FreeMat---" or "---Octave---"
% to find these options and follow the instructions therein.
%%
% Modified 2/28/2018 to include aerodynamic lift
% Ver. 5 changes...
% Modified 3/31/2018 to be compatible with both Octave and Freemat, and
% ...to utilize the odeset Maxstep and Event options.
%%
clear;clc
global Cd Cl0 Clslope rho A S m Pa Pr g
format long
CaseName='';
 while isempty(CaseName)
    CaseName = input('Enter name:','s');
 end
%
V0mph=0;
while (V0mph <= 0)
 V0str='';
 while isempty(V0str)
    V0str = input('Enter V0 in mph:','s'); % Workaround for bug in Freemat V4
 end
V0mph=eval(V0str); % entry speed mph
end
V0=V0mph*0.44704; % entry speed mph m/s
%
theta0deg=0;
while (theta0deg < 5 || theta0deg > 90)
theta0str='';
 while isempty(theta0str)
    theta0str = input('Enter launch degrees (5 to 90):','s'); % Workaround for bug in Freemat V4
 end
theta0deg=eval(theta0str); % initial launch angle in degrees
end
theta0=theta0deg*pi/180; % radians
%
PaDeg=-99;
while (PaDeg < -15 || PaDeg > 15)
PaStr='';
 while isempty(PaStr)
    PaStr = input('Enter attack degrees (-15 to 15):','s'); % Workaround for bug in Freemat V4
 end
PaDeg=eval(PaStr); % Initial attack angle degrees (different than launch angle)
end
Pa=PaDeg*pi/180;
%
PrDeg=-99;
while (PrDeg < -15 || PrDeg > 15)
PrStr='';
 while isempty(PrStr)
    PrStr = input('Enter pitch rate degrees/sec (-15 to 15):','s'); % Workaround for bug in Freemat V4
 end
PrDeg=eval(PrStr); % attack angle rate degrees/sec
end
Pr=PrDeg*pi/180;
%
Y0=2; % Initial height (m)
m=880; % mass of Saab 96, kg (car+passenger 1775lb+165lb)
A=1.795+0.039; % front silhouette area plus extra for hanging wheels, m^2
S=4.146; % surface area of floor, m^2 (approx 44.625 sq ft)
Cd=0.32; % most commonly cited figure
Cl0=0.3; % coefficient of lift at attack angle=0
Clslope=(1.6-Cl0)/(15*pi/180); % coefficient of lift increase per radian (1.6 at 15 degrees)
g=9.81; % m/s^2
rho=1.2041; % density of air, kg/m^3
%% calculate launch velocity based on ramp energy loss of m*g*h
E0=0.5*m*V0^2; % entry energy
Eloss1=m*g*Y0; % ramp height energy loss
%% Further reduce launch energy for drag force acting through length of slope
%% This is a slight overestimate of energy loss because it assumes constant speed
Lslope=Y0/sin(theta0); % Length of slope assuming ramp angle is launch angle
Eloss2=(0.5*Cd*A*rho*V0^2)*Lslope; % energy loss due to drag while on slope
EL=E0-Eloss1-Eloss2; % adjusted launch energy
VL=sqrt(2*EL/m); % launch velocity mps based on launch energy
VLmph=VL/0.44704; % launch velocity mph
%% perform projectile calcs
tmax=5.0; % do calculations for 5 seconds of flight time
tspan=[0 tmax];
%%
% initial conditions as [x0, vx0, y0, vy0]
IC = [0; VL*cos(theta0); Y0; VL*sin(theta0)];
options = odeset('MaxStep',0.01,'Events',@MySaab96JumpEvent);
[t,oput,tf,tfVALS] = ode45(@MySaab96JumpODE, tspan, IC,options);
%%
x= oput(:,1)*3.281; % extract x-position from 1st column
vx= oput(:,2)/0.44704; % extract x-velocity from 2nd column
y= oput(:,3)*3.281; % extract y-position from 3rd column
vy= oput(:,4)/0.44704; % extract y-velocity from 4th column
%%
%%%% Trajectory Calculations Done - Time for Data and Plot Output %%%%
%%
% tf is time at y-->0
xf= tfVALS(1)*3.281; % final x-position at tf
vxf= tfVALS(2)/0.44704; % final x-velocity at tf
yf= tfVALS(3)*3.281; % final y-position at tf
vyf= tfVALS(4)/0.44704; % final y-velocity at tf
%%
Paf=(Pa+Pr*tf)*180/pi % landing attack angle deg
vf = (vxf^2+vyf^2)^0.5 % final velocity
Ymax = max(y) % height of jump at nearest time point
thetaf=atan(vyf/vxf)*180/pi % landing trajectory angle deg
Dy=(sin(30*pi/180)*4+0.833+0.5)*0.3048 % Occupant/driver vertical stopping distance
Gocc=0.5*(vyf*0.44704)^2/Dy/g % (Vyf mph=>m/s)
Dy=(sin(30*pi/180)*4+0.833)*0.3048 % Car vertical stopping distance
Gcar=0.5*(vyf*0.44704)^2/Dy/g % (Vyf mph=>m/s)
Dy=(sin(theta0)*4+0.833+0.5)*0.3048 % Launch perp-to-ramp stopping distance
Vyl=V0*sin(theta0); % velocity perpendicular to ramp m/s
Glaunch=0.5*(Vyl)^2/Dy/g % Launch g on driver
Vylmph=Vyl/0.44704 % velocity perpendicular to ramp mph
%%
figure(1);clf; % default figure options for ---FreeMat--- and ---Octave---
%% To use the following option for ---Octave---, remove the %% and insert %% in the above figure statement...
%% figure(1,'Position',[2,128,1354,320]);clf; % ---Octave--- optional plot resolution [startx,starty,width,height]
%% To use the following option for ---FreeMat---, simply remove the %%...
%% sizefig(1354,320) % ---FreeMat--- optional plot resolution width and height
str=sprintf('%s, V0%6.2f mph, VL%6.2f mph, Angle %4.1f, Attack %4.1f at%4.1f/s, Y0%5.2f ft',CaseName,V0mph,VLmph,theta0deg,PaDeg,PrDeg,Y0*3.281);
if ((PaDeg == 0) && (PrDeg == 0))
  plot(x,y,'b-'); % blue plot if no angle of attack
elseif (PrDeg == 0)
  plot(x,y,'m-'); % magenta plot if constant angle of attack
else
  plot(x,y,'r-'); % red plot if variable angle of attack
end
% NOTE: The following axis statement is optional and is of the form axis([xmin xmax ymin ymax zmin zmax])
% NOTE: Use it only if you want all plots to be of exactly the same scale and axis limits
axis([0 400 0 60 0 1]);
% The following axis equal statement makes sure x and y are spaced equally to give an undistorted plot
axis equal;
grid on;
title(str)
xlabel('Span (ft)')
ylabel('Ht. (ft)')
PlotPng=sprintf('%s.png',CaseName);
print(PlotPng)
FileName=sprintf('%s.txt',CaseName);
str=sprintf('V0:\t\t%6.2f\nVL:\t\t%6.2f\nTheta0:\t\t%6.2f\nPaDeg:\t\t%6.2f\nPrDeg:\t\t%6.2f\ntf:\t\t%6.3f\nxf:\t\t%6.2f\nYmax:\t\t%6.2f\nvxf:\t\t%6.2f\nvyf:\t\t%6.2f\nGlaunch:\t\t%6.2f\nGocc:\t\t%6.2f\nGcar:\t\t%6.2f\nvf:\t\t%6.2f\nPaFinal:\t\t%6.2f\nThetaFinal:\t\t%6.2f\nY0:\t\t%6.2f\n',V0mph,VLmph,theta0deg,PaDeg,PrDeg,tf,xf,Ymax,vxf,vyf,Glaunch,Gocc,Gcar,vf,Paf,thetaf,Y0*3.281);
fp = fopen(FileName,'w');
fprintf(fp,'%s',str);
fclose(fp);
--------------------------------------------------------------------------------------------------------------

File #2, function called by ode45 to update first and second derivatives (velocity and acceleration)
----------------------------------------------- MySaab96JumpODE.m --------------------------------------------
function [ Fout ] = MySaab96JumpODE( t, p )
%% For use with MySaab96Jump.m, Ver. 5, 3/31/2018
%%
%% Copyright 2018, George (Ted) Yurkon
%% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. To
%% view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative Commons,
%% PO Box 1866, Mountain View, CA 94042, USA.
%%
%% In short, this work may be shared provided terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%% International License are adhered to.
%%
% Simultaneous second order differentials for projectile
% motion with air resistance and aerodynamic lift
% modified 2/28/2018 to include aerodynamic lift
% output vector Fout has the four differential outputs x', x'', y', y''
% assumed units: meters, seconds, Newtons, kg, radians
global Cd Cl0 Clslope rho A S m Pa Pr g % these are defined globally so
% they are available to the main script and all functions.
%
% Cdrag=Cd*Afront*rho/2/m drag force constant (Afront varies with attack angle)
% Clift=Cl*Afloor*rho/2/m lift force constant (Cl varies with attack angle)
% 
Pat=Pa+Pr*t; % current attack angle
Afront=A+S*abs(sin(Pat)); % current frontal area includes floor projection due to attack angle
Cdrag=Cd*Afront*rho/2/m;
Cl=Cl0+Clslope*Pat; % current coefficient of lift
Clift=Cl*S*rho/2/m;
Fout = zeros(4,1); % initialize space
Fout(1) = p(2);
Fout(2) = -Cdrag*sqrt(p(2)^2 + p(4)^2)* p(2) - Clift*sqrt(p(2)^2 + p(4)^2)* p(4);
Fout(3) = p(4);
Fout(4) = -g -Cdrag*sqrt(p(2)^2 + p(4)^2)* p(4) + Clift*sqrt(p(2)^2 + p(4)^2)* p(2);
end
--------------------------------------------------------------------------------------------------------------

File #3, an event function which stops ode45 integration when the trajectory crosses the x axis 
----------------------------------------------- MySaab96JumpEvent.m ------------------------------------------
function [value,isterminal,direction] = MySaab96JumpEvent(t,p)% when value is equal to zero, an event is triggered.
%% For use with MySaab96Jump.m, Ver. 5, 3/31/2018
%%
%% Copyright 2018, George (Ted) Yurkon
%% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. To
%% view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative Commons,
%% PO Box 1866, Mountain View, CA 94042, USA.
%%
%% In short, this work may be shared provided terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0
%% International License are adhered to.
%%
% Setting isterminal to 1 stops the solver at the first event.
% Setting isterminal to 0 allows the solver to compute past all events.
% Set direction=0 if all zeros are to be computed (the default)
% Set direction=+1 for only zeros where the event function is increasing.
% Set direction=-1 for only zeros where the event function is decreasing.
value = p(3); % when value = 0 (y=0), an event is triggered
isterminal = 1; % terminate after the first event
direction = -1; % only the zeros when decreasing
end
--------------------------------------------------------------------------------------------------------------

